lambda_seq <- seq(1,.505,-.005)^14
i = 0
lasso.results <- list()
for (year in start.year:end.year) {
  i = i + 1
  print(year)
  # Generate dependent variable
  Violence.mean <-  aggregate(x = as.matrix(dta[dta[,t.var]< year,
                                                     dv]),
                              by = list(dta[dta[,t.var]< year,
                                                 loc.var]),
                              FUN = mean)
  colnames(Violence.mean) <- c(loc.var,
                               "dv_mean")
  Violence.new <- merge(as.matrix(dta[dta[,t.var]< year,
                                           c(loc.var,t.var,dv)]),
                        Violence.mean,
                        by = loc.var)
  Violence.past <- Violence.new[,dv]-Violence.new$dv_mean
  # Fit LASSO
  lasso.results[[i]] <- list()
  dta.past <- dta[dta[,t.var]< year,]
  N <- length(dta.past[,dv])
  cv_results <- foreach (cv=1:cv.runs, 
                         .combine='rbind',
                         .packages = c("glmnet","pROC")) %dopar% {   # Loop over CV runs
                           set.seed(52*cv) 
                           shuffle <- sample(1:N,
                                             N,
                                             replace=FALSE)
                           prediction <- matrix(NA,nrow=N,ncol=length(lambda_seq))
                           sampleSplit = c(0,round((N/5*c(1:5))))
                           for (f in 1:5) {
                             # Fit lasso
                             cv.fit <- glmnet(x = as.matrix(dta.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                                          rhs]),
                                              y = as.matrix(Violence.past[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])]]),
                                              alpha = lasso_demean_params$alpha,
                                              family = lasso_demean_params$family,
                                              lambda = lambda_seq)
                             prediction[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                       1:length(cv.fit$lambda)] <- predict(cv.fit,
                                                                           newx = as.matrix(dta.past[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                                                                          rhs]),
                                                                           type = "response")
                           }
                           cv_results <- c("CV Run"=0,
                                           "lambda"=0,
                                           "MSE"=100)
                           # Get error rates for each lambda used
                           for (l in seq(10,90,10)) {
                             if (sum(is.na(prediction[,l]))==0) {
                               MSE <- mean((Violence.past - prediction[,l])^2)
                               if (MSE<cv_results["MSE"]) {
                                 cv_results <- c("CV Run"=cv,
                                                 "lambda"=lambda_seq[l],
                                                 "MSE" = MSE)
                                 best.l <- l
                               }
                             }
                           }
                           for (l in seq(best.l-9,best.l+10,1)) {
                             if (sum(is.na(prediction[,l]))==0) {
                               MSE <- mean((Violence.past - prediction[,l])^2)
                               if (MSE<cv_results["MSE"]) {
                                 cv_results <- c("CV Run"=cv,
                                                 "lambda"=lambda_seq[l],
                                                 "MSE" = MSE)
                                 best.l <- l
                               }
                             }
                           }
                           return(cv_results)
                         } 
  hours <- floor((proc.time()[3]-start.time)/3600)
  mins <- floor((proc.time()[3]-start.time)/60) - hours*60
  secs <- floor((proc.time()[3]-start.time)) - hours*3600 - mins*60
  
  print(paste(hours,
              "h",
              mins,
              "m",
              secs,
              "s elapsed",
              sep = ""))
  
  optim.lambda <- mean(cv_results[,"lambda"]) 
  lasso.results[[i]]$mod <- glmnet(x = as.matrix(dta[dta[,t.var]< year,
                                                          rhs]),
                                   y = as.matrix(Violence.past),
                                   lambda = seq(10.9,1,-.1)*optim.lambda,
                                   alpha = lasso_demean_params$alpha,
                                   family = lasso_demean_params$family
  )
  opt.lam <- 100
  while (is.na(lasso.results[[i]]$mod$lambda[opt.lam])) {
    opt.lam = opt.lam - 1
  }
  lasso.results[[i]]$lambda <- optim.lambda
  lasso.results[[i]]$fit.oos <- predict(lasso.results[[i]]$mod,
                                        newx = as.matrix(dta[dta[,t.var] == year,
                                                                  rhs]),
                                        type = "response")[,opt.lam]
  Violence.new <- merge(as.matrix(dta[dta[,t.var]== year,
                                           c(loc.var,t.var,dv)]),
                        Violence.mean,
                        by = loc.var)
  lasso.results[[i]]$actual.oos <- Violence.new[,dv]-Violence.new$dv_mean
}
save(lasso.results,file=paste(modeldir,"/",
                              country,
                              "_lasso_",
                              "demean",
                              "_",
                              rhs.group,fileext,
                              ".RData",
                              sep = ""))